Making Maps in R

AnEx Learning Summit
Thursday October 21, 2021

Sarah Rankin
Quantitative Data Analyst
Stategy and Public Impact
New York Public Library
sarahrankin@nypl.org

Topics

• Static maps with sf and ggplot
• Dynamic maps with leaflet

Setup

This is an Rmarkdown notebook. If you’re viewing it in a web browser, you can click on Code->Download Rmd in the top right corner, then open it in RStudio and to run/edit the code. If you don’t use RStudio you can also just copy the code into your editor of choice.

We’ll be using a variety of tidyverse packages as well as sf and leaflet. If you don’t have these packages installed, uncomment the following code and install them now. (It’ll take a little while.)

# install.packages("tidyverse")
# install.packages("sf")
# install.packages("leaflet")

#load just the tidyverse for now
library(tidyverse)

Quick ggplot refresher

Let’s create some very simple point data.

five_points <- data.frame(var1 = 1:5, var2 = as.integer(c(3,5,4,4,1)))
five_points 

• Call ggplot() to make a canvas

• Define aesthetics with aes() - these map your data to spatial properties like the x and y dimensions, or visual aesthetics like color, fill, alpha etc

• Add a geom - here, geom_point()

five_points %>% 
  ggplot(aes(x=var1, y=var2)) +
  geom_point()

• Add flags indicating whether each variable is even or odd (%% is the modulo function)

five_points <- five_points %>%  
  mutate(var1_odd_or_even = ifelse(var1%%2==0,"even","odd"),
         var2_odd_or_even = ifelse(var2%%2==0,"even","odd"))

five_points

• Map “var_1_is even” to the color aesthetic

five_points %>% 
  ggplot(aes(x=var1, y=var2,color = var1_odd_or_even)) +
  geom_point(size = 5) 

Aesthetics in ggplot

• Color is for lines (including outlines of shapes)

• Fill is for areas (filling in shapes)

  • but our points are filled in based on the color aesthetic?!

  • the default point shape is shape 19 (“circle”) - just a solid circle with the color taken from the color aesthetic; shape 21 (“filled circle”) is a circle with a fill and an outline, which can be specified separately

• Generally, you can either map the aesthetic to a data point or specify its value outside of the aes() function. If you do neither, a default will be used.

• Aesthetic mapping/specification can be done in the original ggplot() call, or within individual geoms

Good general reference on specifying aesthetics:

https://ggplot2.tidyverse.org/articles/ggplot2-specs.html#point-1

five_points %>% 
  ggplot(aes(x=var1, 
             y=var2,
             color = var1_odd_or_even,
             fill = var2_odd_or_even)) +
  geom_point(size = 5, 
             shape = 21,
             stroke = 2) 

Points as Maps

• Longitude and Latitude are (on a 2d map) x and y coordinates

• Let’s get some coordinate data from the NYC facilities database

https://data.cityofnewyork.us/City-Government/Facilities-Database/ji82-xba5

public_libraries <- jsonlite::fromJSON("https://data.cityofnewyork.us/resource/ji82-xba5.json?factype=PUBLIC%20LIBRARY") 

public_libraries %>% str()

'data.frame':   221 obs. of  34 variables:
 $ uid       : chr  "015fec883329f3481f1711a39386a5ee" "05375983dd54899915412bc855795f75" "0616b9a31fe2db14793561377be01c0e" "08d71649b166891363a946ebb0b3b029" ...
 $ facname   : chr  "MORRIS PARK LIBRARY" "NORTH FOREST PARK" "58TH STREET LIBRARY" "FLUSHING" ...
 $ addressnum: chr  "985" "98-27" "127" "41-17" ...
 $ streetname: chr  "MORRIS PARK AVENUE" "METROPOLITAN AVENUE" "EAST 58 STREET" "MAIN STREET" ...
 $ address   : chr  "985 MORRIS PARK AVENUE" "98-27 METROPOLITAN AVENUE" "127 EAST 58 STREET" "41-17 MAIN STREET" ...
 $ city      : chr  "BRONX" "FOREST HILLS" "NEW YORK" "FLUSHING" ...
 $ boro      : chr  "BRONX" "QUEENS" "MANHATTAN" "QUEENS" ...
 $ borocode  : chr  "2" "4" "1" "4" ...
 $ zipcode   : chr  "10462" "11375" "10022" "11355" ...
 $ latitude  : chr  "40.84806105240" "40.71118635740" "40.76223213860" "40.75778232920" ...
 $ longitude : chr  "-73.85697658820" "-73.85367190670" "-73.96932374730" "-73.82887651450" ...
 $ xcoord    : chr  "1023819.55661000000" "1024817.14797000000" "992747.99653700000" "1031658.10613000000" ...
 $ ycoord    : chr  "248281.25916200000" "198414.74758100000" "216979.94069500000" "215403.54794900000" ...
 $ bin       : chr  "2045398" "4076687" "1037165" "4114282" ...
 $ bbl       : chr  "2041260040.00000000000" "4032070022.00000000000" "1013130005.00000000000" "4050430011.00000000000" ...
 $ commboard : chr  "211" "406" "105" "407" ...
 $ council   : chr  "13" "29" "4" "20" ...
 $ censtract : chr  "25200" "72900" "11203" "85300" ...
 $ nta       : chr  "BX37" "QN17" "MN19" "QN22" ...
 $ facgroup  : chr  "LIBRARIES" "LIBRARIES" "LIBRARIES" "LIBRARIES" ...
 $ facsubgrp : chr  "PUBLIC LIBRARIES" "PUBLIC LIBRARIES" "PUBLIC LIBRARIES" "PUBLIC LIBRARIES" ...
 $ factype   : chr  "PUBLIC LIBRARY" "PUBLIC LIBRARY" "PUBLIC LIBRARY" "PUBLIC LIBRARY" ...
 $ capacity  : chr  "0" "0" "0" "0" ...
 $ optype    : chr  "Public" "Public" "Public" "Public" ...
 $ opname    : chr  "New York Public Library" "Queens Public Library" "New York Public Library" "Queens Public Library" ...
 $ opabbrev  : chr  "NYPL" "Public" "NYPL" "Public" ...
 $ overlevel : chr  "City" "City" "City" "City" ...
 $ overagency: chr  "New York Public Library" "Queens Public Library" "New York Public Library" "Queens Public Library" ...
 $ overabbrev: chr  "NYPL" "QPL" "NYPL" "QPL" ...
 $ datasource: chr  "nypl_libraries" "qpl_libraries" "nypl_libraries" "qpl_libraries" ...
 $ facdomain : chr  "LIBRARIES AND CULTURAL PROGRAMS" "LIBRARIES AND CULTURAL PROGRAMS" "LIBRARIES AND CULTURAL PROGRAMS" "LIBRARIES AND CULTURAL PROGRAMS" ...
 $ schooldist: chr  "11" "28" "2" "25" ...
 $ policeprct: chr  "49" "112" "18" "109" ...
 $ servarea  : chr  "Local" "Local" "Local" "Local" ...

We need latitude and longitude to be numeric, not strings

public_libraries <- public_libraries %>%  
  mutate(latitude=as.numeric(latitude),longitude=as.numeric(longitude))

public_libraries %>% 
  select(facname,latitude,longitude) %>% 
  str()

'data.frame':   221 obs. of  3 variables:
 $ facname  : chr  "MORRIS PARK LIBRARY" "NORTH FOREST PARK" "58TH STREET LIBRARY" "FLUSHING" ...
 $ latitude : num  40.8 40.7 40.8 40.8 40.7 ...
 $ longitude: num  -73.9 -73.9 -74 -73.8 -73.7 ...

• Map the longitude and latitude to the x and y aesthetics

  • Everyone has a strategy for remembering which is which - I like to think of soup dumplings - “Xiao long bao” - credit to @seankross

public_libraries %>% 
  ggplot(aes(x=longitude,
             y=latitude)) + 
  geom_point() 

• That does not look like a map. Let’s look at our data.

public_libraries %>% 
  select(latitude,longitude) %>% 
  summary()

    latitude       longitude     
 Min.   : 0.00   Min.   :-74.24  
 1st Qu.:40.66   1st Qu.:-73.98  
 Median :40.72   Median :-73.92  
 Mean   :40.35   Mean   :-73.25  
 3rd Qu.:40.77   3rd Qu.:-73.85  
 Max.   :40.90   Max.   :  0.00  

• At least some of our lat/long values are way out of expected range

public_libraries %>% 
  filter(latitude<39|longitude>74) %>% 
  select(latitude,longitude,facname,addressnum,streetname,address,city,boro)

• Looks like we can clean this up by just filtering out these two cases

 public_libraries <- public_libraries %>%   
  filter(!latitude==0,
         !longitude==0)
 
 public_libraries %>% nrow()

[1] 219

• Try again

public_libraries %>% 
  ggplot(aes(x=longitude,
             y=latitude)) + 
  geom_point() 

• Make it a bit prettier

  • Color the points and rename the legend with scale_color_manual()

  • Get rid of axes and background with theme_void()

#define some system colors - a named vector of colors will map the colors to the values of the color aesthetic (backticks let you use non-syntactic R names)
#color values can be names from R's built-in colors or hex codes

library_system_colors <- c(`New York Public Library`="red3",
                   `Brooklyn Public Library` = "orange2",
                   `Queens Public Library` = "purple3")

public_libraries %>% 
  ggplot(aes(x=longitude,
             y=latitude,
             color = overagency)) + 
  geom_point() +
  scale_color_manual(values = library_system_colors,name = "System") +
  theme_void()

The sf (simple features) package

https://r-spatial.github.io/sf/

library(sf)

Simple features is a set of standards for defining two-dimensional geometries, used by various GIS systems, building up from x-y coordinates within a coordinate reference system (crs).

• POINT

• LINESTRING

• POLYGON

• MULTIPOINT

• MULTIPOLYGON

sf package in R lets you:

• store geometry data as a list-column in a data frame or tibble, alongside other data columns

• import data into the R sf format from various GIS formats

• manipulate sf geometries (concatenate, subtract, find overlaps, etc)

You can make an sf object directly from a data frame:

five_points_sf <- five_points %>% st_as_sf(coords = c("var1","var2")) 

five_points_sf %>% print()

Simple feature collection with 5 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1 ymin: 1 xmax: 5 ymax: 5
CRS:           NA
  var1_odd_or_even var2_odd_or_even    geometry
1              odd              odd POINT (1 3)
2             even              odd POINT (2 5)
3              odd             even POINT (3 4)
4             even             even POINT (4 4)
5              odd              odd POINT (5 1)

• There is a special geom for sf objects: geom_sf

• When you add geom_sf to a plot (and the sf package is loaded), ggplot will plot the geometry column - plotting method depends on the sf geometry type (point, polygon, etc)

five_points_sf %>% 
  ggplot() + 
  geom_sf(
    aes(color = var1_odd_or_even,
              fill = var2_odd_or_even),
          size = 5,
          shape = 21,
          stroke = 2
    )

Map Polygons

But for mapping we usually import geometry data from an external dataset.

Formats:

• GEOJSON

• ESRI shapefile

• Many others (st_drivers(what = "vector"))

NYC geography resources:

https://www1.nyc.gov/site/planning/data-maps/open-data/census-download-metadata.page

• Get the NYC PUMA (Public Use Microdata Area) geographies in GEOJSON format (PUMAs in NYC correspond - mostly - to Community Districts)

#read geojson format directly from download link
pumas_2010_geojson <- st_read("https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/NYC_Public_Use_Microdata_Areas_PUMAs_2010/FeatureServer/0/query?where=1=1&outFields=*&outSR=4326&f=pgeojson") %>% 
  st_as_sf()

Reading layer `OGRGeoJSON' from data source 
  `https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/NYC_Public_Use_Microdata_Areas_PUMAs_2010/FeatureServer/0/query?where=1=1&outFields=*&outSR=4326&f=pgeojson' 
  using driver `GeoJSON'
Simple feature collection with 55 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS:  WGS 84

pumas_2010_geojson %>% print()

Simple feature collection with 55 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID PUMA Shape__Area Shape__Length                       geometry
1         1 3701    97928415      53226.65 MULTIPOLYGON (((-73.89641 4...
2         2 3702   188993596     106167.48 MULTIPOLYGON (((-73.86477 4...
3         3 3703   267645154     305260.19 MULTIPOLYGON (((-73.78833 4...
4         4 3704   106217077      47970.15 MULTIPOLYGON (((-73.84793 4...
5         5 3705   122483734      68697.43 MULTIPOLYGON (((-73.87046 4...
6         6 3706    43887042      51825.93 MULTIPOLYGON (((-73.87773 4...
7         7 3707    42281072      37374.60 MULTIPOLYGON (((-73.89964 4...
8         8 3708    55881008      35002.58 MULTIPOLYGON (((-73.92478 4...
9         9 3709   124117798      73287.89 MULTIPOLYGON (((-73.83668 4...
10       10 3710   137760355      90067.74 MULTIPOLYGON (((-73.89681 4...

• Shapefile format is also very common

#for shapefiles, download the zipped shapefile directory into your current working directory, then unzip (can do this outside R if you prefer!)
#download
download.file(url = "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nypuma2010_21c.zip",destfile = "nypuma2010_21c.zip")
#unzip
unzip("nypuma2010_21c.zip")

#read in the shapefile
pumas_2010_shp <- st_read("nypuma2010_21c/nypuma2010.shp") %>% st_as_sf()

Reading layer `nypuma2010' from data source `/Users/sarahrankin/Documents/Mapping/nypuma2010_21c/nypuma2010.shp' using driver `ESRI Shapefile'
Simple feature collection with 55 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 913175.1 ymin: 120121.9 xmax: 1067383 ymax: 272844.3
Projected CRS: NAD83 / New York Long Island (ftUS)

pumas_2010_shp %>% print()

Simple feature collection with 55 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 913175.1 ymin: 120121.9 xmax: 1067383 ymax: 272844.3
Projected CRS: NAD83 / New York Long Island (ftUS)
First 10 features:
   PUMA Shape_Leng Shape_Area                       geometry
1  3701   53227.11   97928517 MULTIPOLYGON (((1012885 268...
2  3702  106167.59  188993632 MULTIPOLYGON (((1021632 267...
3  3703  305269.14  267643637 MULTIPOLYGON (((1042822 243...
4  3704   47970.20  106216909 MULTIPOLYGON (((1026309 256...
5  3705   68697.60  122483690 MULTIPOLYGON (((1020081 255...
6  3706   51826.07   43886931 MULTIPOLYGON (((1018060 261...
7  3707   37374.60   42281072 MULTIPOLYGON (((1012009 253...
8  3708   35002.71   55881068 MULTIPOLYGON (((1005061 247...
9  3709   73288.78  124117768 MULTIPOLYGON (((1029456 237...
10 3710   90067.96  137760290 MULTIPOLYGON (((1012822 229...

• Again, adding geom_sf() to a ggplot renders the data in the geometry column - in this case it draws polygons because the geometry type is multipolygon

pumas_2010_geojson %>% 
  ggplot() + 
  geom_sf()

• Add our points to this map by specifying different data sources within the geom_sf and geom_point calls

• sf objects can have a CRS (coordinate reference system) defined

  ggplot() + 
  geom_sf(data = pumas_2010_geojson,
    fill = "grey",
    color = "grey"
    ) +
  geom_point(data = public_libraries, 
             aes(x=longitude,
                 y=latitude,
                 color = overagency)) +
  scale_color_manual(name = "System",values = library_system_colors) +
  #adding coord sf ensures crs matches
  coord_sf() +
  theme_void()

Chloropleths

• Map areas colored by data value

Broadband adoption in New York City

• NYC Open Data: NYC’s Internet Master Plan: Home Broadband and Mobile Broadband Adoption by PUMA

• Based on American Community Survey Data; shows share of households that have both mobile and home broadband

broadband_use <- read_csv("https://data.cityofnewyork.us/api/views/g5ah-i2sh/rows.csv?accessType=DOWNLOAD")

── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  `PUMA (Public Use Microdata Sample Areas)` = col_double(),
  Borough = col_character(),
  `Home Broadband and Mobile Broadband Adoption (Percentage of  Households)` = col_double(),
  `Home Broadband and Mobile Broadband Adoption by Quartiles (High, Medium-High, Medium-Low, Low)` = col_character()
)

#clean up the names a bit
broadband_use <- broadband_use %>% 
  rename(PUMA = `PUMA (Public Use Microdata Sample Areas)`,
         broadband_adoption = `Home Broadband and Mobile Broadband Adoption (Percentage of  Households)`,
         broadband_adoption_quartile = `Home Broadband and Mobile Broadband Adoption by Quartiles (High, Medium-High, Medium-Low, Low)`)

broadband_use %>% head()

• Add the broadband data to our PUMA geom dataset

#convert PUMA to string; make adoption quartile into factor so it'll sort correctly  
broadband_use <- broadband_use %>% 
  mutate(PUMA = as.character(PUMA),
         broadband_adoption_quartile = factor(broadband_adoption_quartile,
                                              levels = c("High", "Medium High", "Medium Low", "Low")))

pumas_2010_geojson <- pumas_2010_geojson %>% 
  left_join(broadband_use, by = "PUMA")

pumas_2010_geojson %>% print()

Simple feature collection with 55 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID PUMA Shape__Area Shape__Length Borough broadband_adoption broadband_adoption_quartile
1         1 3701    97928415      53226.65   Bronx               0.46                        High
2         2 3702   188993596     106167.48   Bronx               0.30                         Low
3         3 3703   267645154     305260.19   Bronx               0.34                  Medium Low
4         4 3704   106217077      47970.15   Bronx               0.30                         Low
5         5 3705   122483734      68697.43   Bronx               0.32                         Low
6         6 3706    43887042      51825.93   Bronx               0.39                 Medium High
7         7 3707    42281072      37374.60   Bronx               0.42                        High
8         8 3708    55881008      35002.58   Bronx               0.40                 Medium High
9         9 3709   124117798      73287.89   Bronx               0.29                         Low
10       10 3710   137760355      90067.74   Bronx               0.34                  Medium Low
                         geometry
1  MULTIPOLYGON (((-73.89641 4...
2  MULTIPOLYGON (((-73.86477 4...
3  MULTIPOLYGON (((-73.78833 4...
4  MULTIPOLYGON (((-73.84793 4...
5  MULTIPOLYGON (((-73.87046 4...
6  MULTIPOLYGON (((-73.87773 4...
7  MULTIPOLYGON (((-73.89964 4...
8  MULTIPOLYGON (((-73.92478 4...
9  MULTIPOLYGON (((-73.83668 4...
10 MULTIPOLYGON (((-73.89681 4...

• Map the geom_sf fill aesthetic to broadband_adoption

• Use scale_fill_gradient to specify fill colors

  ggplot() + 
  geom_sf(data = pumas_2010_geojson,
    aes(fill = broadband_adoption),
    color = "grey"
    ) +
  geom_point(data = public_libraries,
             aes(x=longitude,
                 y=latitude,
                 color = overagency),
             size = 2) +
  scale_color_manual(name = "System",values = library_system_colors) +
  scale_fill_gradient(low = "grey70",high = "grey10", labels = scales::percent_format(accuracy = 1), name = "Broadband adoption") +
  coord_sf() +
  theme_void()

• Or define discrete colors and use scale_fill_manual

broadband_cols <- c("grey10","grey30","grey50","grey70") %>% set_names(levels(pumas_2010_geojson$broadband_adoption_quartile))


ggplot() + 
  geom_sf(data = pumas_2010_geojson,
    aes(fill = broadband_adoption_quartile),
    color = "grey"
    ) +
  geom_point(data = public_libraries,
             aes(x=longitude,
                 y=latitude,
                 color = overagency),
             size = 2) +
  scale_color_manual(name = "System",values = library_system_colors) +
  scale_fill_manual(values = broadband_cols, name = "Broadband adoption") +
  coord_sf() +
  theme_void()

Median income

https://www1.nyc.gov/site/planning/planning-level/nyc-population/american-community-survey.page

• Download data and add it to our pumas_2010_geojson data frame

download.file(url = "https://www1.nyc.gov/assets/planning/download/office/planning-level/nyc-population/acs/econ_2018_acs5yr_puma.xlsx",
              destfile = "econ_2018_acs5yr_puma.xlsx")

econdata_puma <- readxl::read_xlsx("econ_2018_acs5yr_puma.xlsx", sheet = "EconData")

pumas_2010_geojson <- pumas_2010_geojson %>% 
  left_join(econdata_puma %>% 
              select(PUMA = GeoID,GeogName,median_household_income = MdHHIncE) %>% 
              mutate(median_income_category = 
                       cut(median_household_income,
                           breaks = c(0,50000,75000,100000,150000),
                           labels = c("$0-50K","$50-75K","$75-100K","$100K+"))))

Joining, by = "PUMA"

• Create a new color scale with some greens

median_income_category_cols <- c("#a9ccbc","#7fb29b","#4c7f68","#335545") %>% set_names(levels(pumas_2010_geojson$median_income_category))

ggplot() + 
  geom_sf(data = pumas_2010_geojson,
    aes(fill = median_income_category),
    color = "grey"
    ) +
  geom_point(data = public_libraries,
             aes(x=longitude,
                 y=latitude,
                 color = overagency),
             size = 2) +
  scale_color_manual(name = "System",values = library_system_colors) +
  scale_fill_manual(values = median_income_category_cols, name = "Median income") +
  coord_sf() +
  theme_void() 

Leaflet: interactive maps

library(leaflet)

Leaflet: open source javascript library for making interactive maps on various platforms

R has a leaflet package that lets you create leaflet “widgets” within R

Create a map, add tiles, set the view

• Similar to ggplot in that we keep adding layers/elements, except we continue using the magrittr pipe %>% to add elements to a leaflet map (instead of the ggplot +)

• for now we’re creating a map without any data

leaflet() %>% 
  #addTiles() %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lat = 40.7,lng = -74, zoom = 10)

Add our public library branches

• addCircles() adds points to the map

• can specify the data source within the addCircles() call or in the original leaflet() call

• to refer to fields in our data frame, use the tilde prefix ~

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lat = 40.7,lng = -74, zoom = 10) %>%
  addCircles(data = public_libraries,
             lng = ~longitude,
             lat = ~latitude
             )

• Use colorFactor() to create a palette to color the circles by system

• Use the label and popup arguments to add hover/click info

pal_system <- colorFactor(palette = library_system_colors ,
                          domain = unique(public_libraries$overagency),
                          ordered = T)

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lat = 40.7,lng = -74, zoom = 10) %>%
  addCircles(data = public_libraries,
             lng = ~longitude,
             lat = ~latitude,
             color = ~pal_system(overagency),
             fill = ~pal_system(overagency),
             radius = 100,
             label = ~facname,
             popup = ~address
             ) %>% 
  addLegend(position = "topleft",
            pal = pal_system, 
            values = unique(public_libraries$overagency))

Because we have a basemap, we don’t need the PUMA polygons to create a basic map. But we can add them if we want our chloropleths on the interactive map.

• Fill color argument in leaflet is fillColor; specify opacity with fillOpacity

• Use stroke to specify whether boundary lines are drawn; color , weight, and opacity to control their appearance

pal_broadband <- colorFactor(palette = broadband_cols, 
                             domain = levels(pumas_2010_geojson$broadband_adoption_quartile),
                             ordered = T)


leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lat = 40.7,lng = -74, zoom = 10) %>%
  addPolygons(data = pumas_2010_geojson,
              fillColor = ~pal_broadband(broadband_adoption_quartile),
              fillOpacity = .7,
              color = "grey",
              stroke = T,
              weight = 1,
              label = ~PUMA,
              popup = ~paste0("<b>PUMA ",PUMA,"</b>",
                             "<br>", GeogName,
                             "<br>Broadband adoption: ",round(broadband_adoption*100,0),"%",
                             "<br>Broadband adoption category: ",broadband_adoption_quartile)) 

Combine it with the libraries and add a legend

• Order matters - each successive layer is plotted on top of the previous ones

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lat = 40.7,lng = -74, zoom = 10) %>%
  addPolygons(data = pumas_2010_geojson,
              fillColor = ~pal_broadband(broadband_adoption_quartile),
              fillOpacity = .7,
              color = "grey",
              stroke = T,
              weight = 1,
              label = ~PUMA,
              popup = ~paste0("<b>PUMA ",PUMA,"</b>",
                              "<br>", GeogName,
                             "<br>Broadband adoption: ",round(broadband_adoption*100,0),"%",
                             "<br>Broadband adoption category: ",broadband_adoption_quartile)) %>% 
    addCircles(data = public_libraries,
             lng = ~longitude,
             lat = ~latitude,
             color = ~pal_system(overagency),
             fill = ~pal_system(overagency),
             radius = 100,
             label = ~facname,
             popup = ~address
             ) %>% 
  addLegend(position = "topleft",
            pal = pal_system, 
            values = public_libraries$overagency,
            title = "Library System") %>% 
  addLegend(position = "topleft",
            pal = pal_broadband, 
            values =  pumas_2010_geojson$broadband_adoption_quartile,
            opacity = .7,
            title = "Broadband Access")

Visualize more than one layer

• Use group argument within each layer to specify which group it belongs to

• addLayersControl adds a control to toggle between layers or turn layers on and off

• assigning the map to an object creates a leaflet object with class htmlwidget

pal_income <- colorFactor(palette = median_income_category_cols, 
                             domain = levels(pumas_2010_geojson$median_income_category),
                             ordered = T)

libraries_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lat = 40.7,lng = -74, zoom = 10) %>%
  addPolygons(data = pumas_2010_geojson,
              fillColor = ~pal_broadband(broadband_adoption_quartile),
              fillOpacity = .7,
              color = "grey",
              stroke = T,
              weight = 1,
              label = ~PUMA,
              popup = ~paste0("<b>PUMA ",PUMA,"</b>",
                              "<br>", GeogName,
                             "<br>Broadband adoption: ",round(broadband_adoption*100,0),"%",
                             "<br>Broadband adoption category: ",broadband_adoption_quartile),
              group = "Broadband") %>% 
  addPolygons(data = pumas_2010_geojson,
              fillColor = ~pal_income(median_income_category),
              fillOpacity = .7,
              color = "grey",
              stroke = T,
              weight = 1,
              label = ~PUMA,
              popup = ~paste0("<b>PUMA ",PUMA,"</b>",
                              "<br>", GeogName,
                             "<br>Median income: $",format(median_household_income,big.mark = ","),
                             "<br>Median income category: ",median_income_category),
              group = "income") %>% 
    addCircles(data = public_libraries,
             lng = ~longitude,
             lat = ~latitude,
             color = ~pal_system(overagency),
             fill = ~pal_system(overagency),
             opacity = 1,
             fillOpacity = 1,
             radius = 100,
             label = ~facname,
             popup = ~address#,
             #group = "Broadband"
             ) %>% 
    # addCircles(data = public_libraries,
    #          lng = ~longitude,
    #          lat = ~latitude,
    #          color = ~pal_system(overagency),
    #          fill = ~pal_system(overagency),
    #          opacity = 1,
    #          fillOpacity = 1,
    #          radius = 100,
    #          label = ~facname,
    #          popup = ~address,
    #          group = "Income"
    #          ) %>% 
  addLegend(position = "topleft",
            pal = pal_system, 
            values = public_libraries$overagency,
            title = "Library System") %>% 
  addLegend(position = "topleft",
            pal = pal_broadband, 
            values =  pumas_2010_geojson$broadband_adoption_quartile,
            opacity = .7,
            title = "Broadband Access",
            group = "Broadband") %>% 
  addLegend(position = "topleft",
            pal = pal_income, 
            values =  pumas_2010_geojson$median_income_category,
            opacity = .7,
            title = "Median household income",
            group = "Income") %>% 
  addLayersControl(
            baseGroups = c("Broadband","Income"),
            options = layersControlOptions(collapsed = F)#,
            #overlayGroups = c("Branches")
  ) 


libraries_map

The htmlwidgets package lets you save this map as an html file, which can then be opened in a browser or embedded in a web page

htmlwidgets::saveWidget(libraries_map,
                        file = "libraries_map.html",
                        title = "NYC Libraries:  Broadband Access and Median Income")

Other resources

• tidycensus package: amazing R interface for Census/ACS API, including geometries

  • Analyzing Census Data: Methods Maps and Models in R

• tigris package: more Census geometries

• rnaturalearth: world map data

• tmap: Alternative package for tmapping in R, with ggplot-like syntax - some find this the best way to get started

• mapboxapi package: more basemap options, geocoding services, isochrones (requires setting up a mapbox account)

• Geocoding

  • Within NYC, Geosupport

  • Urban Institute overview of other tools

• Geocomputation with R